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We revisit the classical stability versus accuracy dilemma for the lattice Boltzmann methods 
(LBM). Our goal is a stable method of second-order accuracy for fluid dynamics based on the 
lattice Bhatnager-Gross-Krook method (LBGK). 

The LBGK scheme can be recognised as a discrete dynamical system generated by free-flight and 
entropic involution. In this framework the stability and accuracy analysis are more natural. We 
find the necessary and sufficient conditions for second-order accurate fluid dynamics modelling. In 
particular, it is proven that in order to guarantee second-order accuracy the distribution should 
belong to a distinguished surface - the invariant film (up to second-order in the time step). This 
surface is the trajectory of the (quasi) equilibrium distribution surface under free-flight. 

The main instability mechanisms are identified. The simplest recipes for stabilisation add no 
artificial dissipation (up to second-order) and provide second-order accuracy of the method. Two 
other prescriptions add some artificial dissipation locally and prevent the system from loss of posi- 
tivity and local blow-up. Demonstration of the proposed stable LBGK schemes are provided by the 
numerical simulation of a ID shock tube and the unsteady 2D-flow around a square-cylinder up to 
Reynolds number 0(10000). 



I. INTRODUCTION 

A lattice Boltzmann method (LBM) is a discrete ve- 
locity method in which a fluid is described by associat- 
ing, with each velocity Vi, a single-particle distribution 
function /j = fi(x,t) which is evolved by advection and 
interaction on a fixed computational lattice. 

The method has been proposed as a discretization of 
Boltzmann's kinetic equation (for an introduction and 
historic review see [53| ). Furthermore, the collision oper- 
ator can be alluringly simplified, as is the case with the 
Bhatnager-Gross-Krook (BGK) operator Q, whereby 
collisions are described by a single-time relaxation to lo- 
cal equilibria /*: 

df 1 

Jl +Vi . Vfi = -(/*-/,). (1) 

The physically reasonable choice for /* is as entropy max- 
imizers, although other choices of equilibria are often 
preferred [53|. The local equilibria /* depend nonlin- 
early on the hydrodynamic moments (density, momen- 
tum, etc.). These moments are linear functions of /i, 
hence ([TJ is a nonlinear equation. For small r, the 
Chapman-Enskog approximation [l4[ reduces (|T|) to the 
compressible Navier-Stokes equation [53| with kinematic 
viscosity v ~ tc\, where c\ is the thermal velocity for 
one degree of freedom. 

The overrelaxation discretization of (P) (see, e.g., [f| 
Gil US US SO, Hl| ) is known as LBGK, and allows one to 
choose a time step St r. This decouples viscosity from 
the time step, thereby suggesting that LBGK is capa- 
ble of operating at arbitrarily high- Reynolds number by 
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making the relaxation time sufficiently small. However, 
in this low-viscosity regime, LBGK suffers from numeri- 
cal instabilities which readily manifest themselves as local 
blow-ups and spurious oscillations. 

Another problem is the degree of accuracy. An approx- 
imation to the continuous-in-time kinetics is not equiva- 
lent to an approximation of the macroscopic transport 
equation. The fluid dynamics appears as a singular 
limit of the Boltzmann or BGK equation for small r. 
An approximation to the corresponding slow manifold in 
the distribution space is constructed by the Chapman- 
Enskog expansion. This is an asymptotic expansion, and 
higher (Burnett) terms could have singularities. An alter- 
native approach to asymptotic expansion (with "diffusive 
scaling" instead of "convective scaling" in the Chapman- 
Enskog expansion) was developed in |37j in order to ob- 
tain the incompressible Navier-Stokes equations directly 
from kinetics. 

It appears that the relaxation time of the overrelax- 
ation scheme to the slow hydrodynamic manifold may 
be quite large for small viscosity: t rc i ax ~ c\5f/(2v) ~ 
5t/(2r) (see below, in Sec. Some estimates of long 
relaxation time for LBGK at large Reynolds number are 
found earlier in [58[. So, instead of fast relaxation to 
a slow manifold in continuous-in-time kinetics, we could 
meet a slow relaxation to a fluid dynamics manifold in 
the chain of discrete LBM steps. 

Our approach is based on two ideas: the Ehrenfests' 
coarse-graining [13, HH, [3(| and the method of differen- 
tial approximation of difference equations [34], [5fjj . The 
background knowledge necessary to discuss the LBM in 
this manner is presented in Sect. [TTJ In this section, we 
answer the question: how to provide second-order accu- 
racy of the LBM methods for fluid dynamics modelling? 
We prove the necessary and sufficient conditions for this 
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accuracy. It requires a special connection between the 
distribution fa and the hydrodynamic variables. There 
is only one degree of freedom for the choice of fi, if the 
hydrodynamic fields are given. Moreover, the LBM with 
overrelaxation can provide approximation of the macro- 
scopic equation even when it does not approximate the 
continuous-in-time microscopic kinetics. 

This approach suggests several sources of numerical 
instabilities in the LBM and allows several recipes for 
stabilisation. A geometric background for this analysis 
provides a manifold that is a trajectory q of the quasiequi- 
librium manifold due to free-flight. We call this manifold 
the invariant film (of nonequilibrium states). It was in- 
troduced in and studied further in [H, H3, HI] . Com- 
mon to each stabilisation recipe is the desire to stay uni- 
formly close to the aforementioned manifold (Sect. IIII[) . 

In Sect. HV| in addition to two LBM accuracy tests, 
a numerical simulation of a ID shock tube and the 
unsteady 2D-flow around a square-cylinder using the 
present stabilised LBM are presented. For the later prob- 
lem, the simulation quantitatively validates the experi- 
mentally obtained Strouhal-Reynolds relationship up to 
Re = 0(10000). This extends previous LBM studies of 
this problem where the relationship had only been suc- 
cessfully validated up to Re = 0(1000) 

Sect. [V] contains some concluding remarks as well as 
practical recommendations for LBM realisations. 

We use operator notation that allows us to present gen- 
eral results in compact form. The only definition we have 
to recall here is the (Gateaux) differential: the differential 
of a map </(/) at a point /o is a linear operator (DfJ) f 
defined by a rule: (D f J) fo g = ^(J(/o + eg)) e=0 . 

II. BACKGROUND 



state, and M for the macroscopic state. The vector M is 
a linear function of /: M — m(f). 

b. Equilibrium. For any allowable value of M an 
"equilibrium" distribution should be given: a microscopic 
state flj. It should satisfy the obvious, but important 
identity of self-consistency: 

m(f* M ) = M, (2) 

or in differential form 

m(D M f*M) = 1, i-e-, m((D M flj)a) = a. (3) 

The state is not a proper thermodynamic equilibrium, 
but a conditional one under the constraint m(f) = M. 
Therefore we call it a quasiequilibrium (other names, such 
as local equilibrium, conditional equilibrium, generalised 
canonical state or pseudoequilibrium are also in use). 

For the quasiequilibrium an equilibration opera- 
tion is the projection II* of the distribution / into the 
corresponding quasiequilibrium state: II*(/) = f^tfy 

In the fully physical situation with continuous veloc- 
ity space, the quasiequilibrium is defined as a condi- 
tional entropy maximizer by a solution of the optimisa- 
tion problem: 

S(J) - max, m(f) = M, (4) 

where S(f) is an entropy functional. 

The choice of entropy is ambiguous; generally, we can 
start from a concave functional of the form 

S(f) = f s(f(x, v, t))f(x, v, t) dxdv (5) 



a. Microscopic and macroscopic variables. Let us 
describe the main elements of the LBM construction. 
The first element is a microscopic description, a single- 
particle distribution function f(x,v), where x is the 
space vector, and v is velocity. If velocity space is approx- 
imated by a finite set {vi}, then the distribution is ap- 
proximated by a measure with finite support, f(x, v) sw 
J2i fiS( v — v i)- In that case, the microscopic description 
is the finite-dimensional vector- function fi(x). 

The second main element is the macroscopic descrip- 
tion. This is a set of macroscopic vector fields that 
are usually some moments of the distribution func- 
tion. The main example gives the hydrodynamic fields 
(density-momcntum-energy density): {n,nu, £}(x) = 
/{l, v, v 2 /2}f(x, v) dv. But this is not an obligatory 
choice. If we would like to solve by LBM methods the 
Grad equations l3ll. |49| or some extended thermody- 
namic equations [361 ] . we should extend the list of mo- 
ments (but, at the same time, we should be ready to in- 
troduce more discrete velocities for a proper description 
of these extended moment systems). 

In general, we use the notation / for the microscopic 



with a concave function of one variable s(f). The choice 
by default is s(f) = —In/, which gives the classical 
Boltzmann-Gibbs-Shannon (BGS) entropy. 

For discrete velocity space, there exist some extra mo- 
ment conditions on the equilibrium construction: in ad- 
dition to (J2) some higher moments of a discrete equilib- 
rium should be the same as for the continuous one. This 
is necessary to provide the proper macroscopic equations 
for M. Existence of entropy for the entropic equilibrium 
definition (01 whilst fulfilling higher moment conditions 
could be in contradiction, and a special choice of velocity 
set may be necessary (for a very recent example of such 
research for multispeed lattices see [l?} )• Another choice 
is to refuse to deal with the entropic definition of equi- 
librium (|4|) and assume that there will be no perpetuum 
mobile of the second kind. This extends the possibility 
for approximation, but creates some risk of nonphysical 
behavior of the model. For a detailed discussion of the 
-H-theorem for LBM we refer the readers to (54|. 

Some of the following results depend on the entropic 
definition of equilibrium, but some do not. We always 
point out if results are "entropy-free" . 
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c. Free flight. In the LBM construction the other 
main elements are: the free-flight transformation and the 
collision. There are many models of collisions, but the 
free-flight equation is always the same 



df 



(0) 



with exact solution f(x,v,t) = f(x — vt,v,0), or for 
discrete velocities, 



dfi 
dt 



+ Vi ■ VJ, = 0, 



(7) 



fi(x, t) = fi(x — Vit, 0). Free-flight conserves any entropy 
of the form |5]). In general, we can start from any dy- 
namics. For application of the entropic formalism, this 
dynamics should conserve entropy. Let this kinetic equa- 
tion be 



(8) 



For our considerations, the free-flight equation will be the 
main example of the conservative kinetics (j8]) . 

The phase flow 9 t for kinetic equation © is a shift in 
time that transforms /(to) m to /(to + t). For free-flight, 
6 t : f(x,v) -> f(x - vt,v). 

Remark. We work with dynamical systems defined 
by partial differential equations. Strictly speaking, this 
means that the proper boundary conditions are fixed. 
In order to separate the discussion of equation from a 
boundary condition problem, let us imagine a system 
with periodic boundary conditions (e.g., on a torus), or a 
system with equilibrium boundary conditions at infinity. 

d. Ehrenfests ' solver of second-order accuracy for the 
Navier-Stokes equations. Here we present a generali- 
sation of a well known result. Let us study the fol- 
lowing process (an example of the Ehrenfests' chain 
[22l 124 . [30( | . a similar result gives the optimal predic- 
tion approach (l8j): free-flight for time r - equilibration 
- free-flight for time r - equilibration - ■ ■ ■ . During this 
process, the hydrodynamic fields approximate the solu- 
tion of the (compressible) Navier-Stokes equation with 



viscosity v 
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where c\ is the thermal velocity for 



one degree of freedom. The error of one step of this ap- 
proximation has the order C(r 3 ). An exact expression 
for the transport equation that is approximated by this 
process in the general situation (for arbitrary initial ki- 
netics, velocity set and for any set of moments) is: 

AM t 

— = m{J c {f M )) + -m((D f J c (f))f L A r J, (9) 

where A f* is the defect of invariance of the quasiequi- 
librium manifold: 

A /Jr = j c (/m) - D M (r M )m(j c (r M )), (10) 

and is the difference between the vector-field J c and its 
projection on to the quasiequilibrium manifold. This re- 
sult is entropy-free. 



The first term in the right hand side of ([9]) - the 
quasiequilibrium approximation - consists of moments of 
df/dt computed at the quasiequilibrium point. For free- 
flight, hydrodynamic fields and Maxwell equilibria this 
term gives the Euler equations. The second term in- 
cludes the action of the differential DfJ c (f)f* on the 
defect of invariance A/« (for free-flight ©, this differ- 
ential is just —v ■ W x , for the discrete version (JT]) this 
is the vector-column —v. L ■ V x ). These terms always ap- 
pear in the Chapman-Enskog expansion. For free-flight, 
hydrodynamic fields and Maxwell equilibria they give 
the Navier-Stokes equations for a monatomic gas with 
Prandtl number Pr = 1: 
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where m is particle mass, ks is Boltzmann's constant, 
P = nk^T is ideal gas pressure, T is kinetic temperature, 
and the underlined terms are the results of the coarse- 
graining additional to the quasiequilibrium approxima- 
tion. 

All computations are straightforward exercises (differ- 
ential calculus and Gaussian integrals for computation of 
the moments, m, in the continuous case). More details 
of these computations are presented in [28[ . 

The dynamic viscosity in (fTTj) is fx = ^nk^T. It is use- 
ful to compare this formula to the mean-free-path the- 
ory that gives /i = r co \nk^,T, where t co \ is the collision 
time (the time for the mean- free-path) . According to 
these formulas, we get the following interpretation of the 
coarse-graining time r for this example: r = 2r co i. 

For any particular choice of discrete velocity set {v^ 
and of equilibrium f^ the calculation could give different 
equations, but the general formula ([9]) remains the same. 
The connection between discretization and viscosity was 
also studied in [5l| . Let us prove the general formula 

We are looking for a macroscopic system that is ap- 
proximated by the Ehrenfests' chain. Let us look for 
macroscopic equations of the form 



dM 
~dT 



*(M) 



(12) 



with the phase flow 4> t : M(t) = $ t M(0). The trans- 
formation $ T should coincide with the transformation 
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M i ^ m(0 r (/L)) up to second-order in r. The match- 
ing condition is 

m(e T (/; / )) = $ r (M) for every M and given r. (13) 

This condition is the equation for the macroscopic vector 
field *if(M). The solution of this equation is a function of 
r: \J> — \I>(M, r). For a sufficiently smooth microscopic 
vector field J c (f) and entropy it is easy to find the 
Taylor expansion of ^(M, r) in powers of r. Let us find 
the first two terms: *(M,r) = * (M) +r*i(M) +o(r). 
Up to second-order in r the matching condition ()13|) is 



rn(J c (f* M ))T + m((D f J c (f)) f * M (J c (r M )))- 



= * (M)r + *i(M)r 2 + (J) w $ (itf))($ (M))- 



From this condition immediately follows: 



(14) 



* (M) = m{Ur M % 

*i(M) = im((D/J c (/)) /5f A /Jf ), 



where A ^» is the defect of invariance (|T0|) . Thus we find 
that the macroscopic equation in the first approximation 
is ©. 

e. The Chapman-Enskog expansion for the gener- 
alised BGK equation. Here we present the Chapman- 
Enskog method for a class of generalised model equa- 
tions. This class includes the well-known BGK kinetic 
equation, as well as many other model equations [26| . 

As a starting point we take a formal kinetic equation 
with a small parameter r 

¥ = Af)~ J c (/) + -(n*(/)-/). (i6) 

at t 

The term IT(/) — / is nonlinear because of the nonlinear 
dependency of IT(/) = f* n{f) on m(f). 

We would like to find a reduced description valid for 
the macroscopic variables M. This means, at least, that 
we are looking for an invariant manifold parameterised 
by M, f — /m, that satisfies the invariance equation: 



(D M fM)(m(J(f M ))) = J(fu). 



(17) 



The invariance equation means that the time deriva- 
tive of / calculated through the time derivative of M 
(M = to(J(/m))) by the chain rule coincides with the 
true time derivative J(/). This is the central equation 
for model reduction theory and applications. The first 
general results about existence and regularity of solutions 
to (fT7| were obtained by Lyapunov [45| (see, e.g., the re- 
view in [HI]). For the kinetic equation (|16p the invariance 
equation has the form 



(D M f M )(m(J c (f M ))) = J c (/m) + -(f* M - f M ), (18) 

r 



because of the self-consistency identity d2j) , (J3J) . 

Due to the presence of the small parameter t in J(f), 
the zeroth approximation to /m is the quasiequilibrium 

f(0) 



approximation: f A[ = f^- Let us look for /m in the 
form of a power series: /m = /£ + t /m'' + ' ' ' j with 
m(f^ ) = for k > 1. From (fTS]) we immediately find: 



/m = U$) - {D M f^){m{J c (f^))) = A /if . (19) 

It is very natural that the first term of the Chapman- 
Enskog expansion for the model equation (fTo]) is just the 
defect of invariance for the quasiequilibrium (jlOp . 

The corresponding first-order in r approximation for 
the macroscopic equations is: 

^ = m{J c {f* M )) + rm((D f J c U))flAfh)- ( 20 ) 



We should recall that m(A/j f ) = 0. The last term in l|18]) 
vanishes in the macroscopic projection for all orders. The 
only difference between and © is the coefficient 1/2 
before r in ©. 

/. Decoupling of time step and viscosity: how to pro- 
vide second-order accuracy? In the Ehrenfests' chain 
"free-flight - equilibration - • • • " the starting point of 
each link is a quasiequilibrium state: the chain starts 
from /^f( ), then, after free-flight, equilibrates into f^^y 
etc. The viscosity coefficient in is proportional to r. 
Let us choose another starting point f^ f in order to de- 
couple time step and viscosity and preserve the second- 
order accuracy of approximation. We would like to get 
equation ([9]) with a chain time step St — h. Analogously 
to Q14j) and (fl~5)) , we obtain the macroscopic equation 

^ - m(J c (f M ))+m((D f J c (f)) rM ((/WM)+(V2)A /Jtf )), 

(21) 

under the condition that f^ — f^ = O(h). The initial 
point 



f'li — f. 



M 



r)^r M + o{h) 



(22) 



provides the required viscosity. This is a sufficient con- 
dition for the second-order accuracy of the approxima- 
tion. Of course, the self-consistency identity m(flj) = M 
should be valid exactly, as (J2J) is. This starting distribu- 
tion is a linear combination of the quasiequilibrium state 
and the first Chapman-Enskog approximation. 

The necessary and sufficient condition for second-order 
accuracy of the approximation is: 

m ((D f J c (f)) r jr M - f* M + \{h r)A /Jf )) = o(h) 

(23) 

(with the self-consistency identity m(/|, f ) = M). This 
means, that the difference between left and right hand 
sides of (|2"2")l should have zero moments and give zero 
inputs in observable macroscopic fluxes. 
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Hence, the condition of second-order accuracy signif- 
icantly restricts the possible initial point for free- flight. 
This result is also entropy-free. 

Any construction of collisions should keep the system's 
starting free-flight steps near the points given by ([22]) 
and (|23|) . The conditions (|22|) and (f23|) for second-order 
accuracy of the transport equation approximation do not 
depend on a specific collision model, they are valid for 
most modifications of the LBM and LBGK that use free- 
flight as a main step. 

Various multistep approximations give more freedom 
of choice for the initial state. For the construction of 
such approximations below, the following mean viscosity 
lemma is important: if the transformations f2|* : M — > 
M', i = 1, . . . , k, approximate the phase flow for © for 
time h (shift in time h) and r = n with second-order ac- 
curacy in h, then the superposition f^Slj • • • approx- 
imates the phase flow for (© for time kh (shift in time 
kh) for the average viscosity r = x(ti + ■ ■ ' + 17.) with the 
same order of accuracy. The proof is by straightforward 
multiple applications of Taylor's formula. 

g. Entropic formula for D m (/X/ ) ■ Among the many 
benefits of thermodynamics for stability analysis there 
are some technical issues too. The differential of equilib- 
rium £*m(/m) appears in many expressions, for example 
©, CLU]), dill), (IS]), (HZ]) and (TTH]). If the quasiequihb- 
rium is defined by the solution of the optimisation prob- 
lem (gj, then 

D M f* M = (DjSy^ m T (m (DjS)'^ m T ) ~* . (24) 

This operator is constructed from the vector m, the 
transposed vector m T and the second differential of en- 
tropy. The inverse Hessian (d 2 S/dftdfj)^ 1 is especially 
simple for the BGS entropy, it is just foSij. The for- 
mula ([24]) was first obtained in [48| (for an important 
particular case; for further references see [28|). 

h. Invariant film. All the points Sj(/L-) belong to 
a manifold that is a trajectory q of the quasiequihbrium 
manifold due to the conservative dynamics (JSJ) (in hydro- 
dynamic applications this is the free-flight dynamics l[5]). 
We call this manifold the invariant film (of nonequilib- 
rium states). It was introduced in [25| and studied fur- 
ther in [HEzLSl]. The defect of invariance A/j, (TDJ) is 
tangent to q at the point /j^, and belongs to the intersec- 
tion of this tangent space with kerm. This intersection is 
one-dimensional. This means that the direction of Ap 
is selected from the tangent space to q by the condition: 
derivative of M in this direction is zero. 

A point / on the invariant film q is naturally param- 
eterised by (M, t): f = qu,t, where M = m(f) is the 
value of the macroscopic variables, and t — t(f) is the 
time shift from a quasiequihbrium state: 0_ t (/) is a 
quasiequihbrium state for some (other) value of M. By 
definition, the action of @t on the second coordinate of 
qm.t is simple: Qt(qM,r) = QM',t+r- To the first-order in 
t, ' 



and (?m,o = /m- The quasiequihbrium manifold divides 
q into two parts, q = q_ U qo U q+, where q_ = {qu,t\ t < 
0}, q+ = {qM,t\t > 0}, and qo is the quasiequihbrium 
manifold: q = Um,o} = {fti}- 

There is an important temporal involution of the film: 

Ir{<lM,t) = QM,-t- (26) 

Due to (|22[) . for qu,t and a given time step h the trans- 
formation M <— ► m(Oh{qM,t)) approximates the solution 
of with t = It + h for the initial conditions M and 
time step h with second-order accuracy in h. Hence, due 
to the mean viscosity lemma, the two-step transforma- 
tion 

M i * m(I T (e h (I T (e h (q M ,t))))) (27) 

approximates the solution of §§§ with r = (the Eulcr 
equations) for the initial conditions M and time step h 
with second-order accuracy in h. This is true for any t, 
hence, for any starting point on the invariant film with 
the given value of M. 

To approximate the solution of © with nonzero r, we 
need an incomplete involution: 

IxilMt) = 9M,-(2/3-l)t- (28) 

For j3 = 1, we have l\ — It and for (3 — 1/2, I^/ 2 is 
just the projection onto the quasiequihbrium manifold: 

1 /2 

Iji (qM,t) = n*(<jM,t) = <Zm,o- After some initial steps, 
the following sequence gives a second-order in time step h 
approximation of © with r = (1 - (3)h/(3, 1/2 < /3 < 1: 

M n = m((4e h ) n q M ,t)- (29) 

To prove this statement we consider a transformation of 
the second coordinate in qM,s n by I^Qh'- 

d n+1 = -{2fi-l)(d n + h). (30) 

This transformation has a fixed point $* = —h(2p — 
l)/(2/3) and i? n = + {-l) n {2(3 - l) n 5 for some S. If 
1 — f3 is small then relaxation may be very slow: d n ~ 
■d* + (— l) n 5 exp(— 2n(l — /3)), and relaxation requires ~ 
1/(2(1 — P)) steps. If i!} n = ??* + o(h) then the sequence 
M n (J21D approximates © with r = h-2\&*\ = (l-/3)ft//3 
and second-order accuracy in the time step h. The fixed 
points (?a/,iS* coincide with the restart points f^j+^Af^ 
(|22| in the first order in d* — —(h — r)/2, and the middle 
points i9* +h/2 of the free-flight jumps <7a/,i?* i_ ► qM',$*+h 
approximate the first-order Chapman-Enskog manifold 

For the entropic description of quasiequihbrium, we 
can connect time with entropy and introduce entropic 
coordinates. For each M and positive s from some in- 
terval < s < c there exist two numbers t± (M, s) 
(f+(M, s) > 0, i_(M, s) < 0) such that 



(25) 



S(qM,t±(M,s)) = S{/m) - s. (31) 
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The numbers t± coincide to the first-order: t + = — t_ + 
(t_). 

We define the entropic involution Is as a transforma- 
tion of q: 

is(m,t±) = iM,t T - (32) 

The introduction of incomplete entropic involution Ig is 
also obvious (see [241]). 

Entropic involution Ig coincides with the temporal in- 
volution It, up to second-order in the deviation from 
quasiequilibrium state / — n*(/). Hence, in the vicin- 
ity of quasiequilibrium there is no significant difference 
between these operations, and all statements about the 
temporal involution are valid for the entropic involution 
with the same level of accuracy. 

For the transfer from free-flight with temporal or en- 
tropic involution to the standard LBGK models we must 
transfer from dynamics and involution on q to the whole 
space of states. Instead of Ij, or I<t the transformation 

4 :/_► H*(/) + (2/? -1)(II* (/)-/) (33) 
is used. For f3 — 1, Iq is a mirror reflection in the 

1 /2 

quasiequilibrium state IT*(/), and for j3 — 1/2, I ' is 
the projection onto the quasiequilibrium manifold. If, for 
a given f = qu,t, the sequence (|2Uj) gives a second-order 
in time step h approximation of ([9]), then the sequence 

M n = m((I^e h ) n fo) (34) 

also gives a second-order approximation to the same 
equation with r = (1— f3)h/f3 . This chain is the standard 
LBGK model. 

Entropic LBGK (ELBGK) methods d HIE 13 dif- 
fer only in the definition of (f3"3"|) : for (3 = 1 it should 
conserve entropy, and in general has the following form: 

/|(/) = (l-/3)/ + Af, (35) 

with / = (1 — a)f + all*(f). The number a = a(f) is 
chosen so that a constant entropy condition is satisfied: 
S(f) = S(f). For LBGK (J33j>, a = 2. 

Of course, computation of Iq is much easier than that 

of ly, I a or A: it is not necessary to follow exactly the 
manifold q and to solve the nonlinear constant entropy 
condition equation. For an appropriate initial condition 
from q (not sufficiently close to qo), two steps of LBGK 
with Iq give the same second-order accuracy as ([2T)| . But 
a long chain of such steps can lead far from the quasiequi- 
librium manifold and even from q. Here, we see stability 
problems arising. For (3 close to 1, the one-step transfor- 
mation IqQii in the chain ([34]) almost conserves distances 
between microscopic distributions, hence, we cannot ex- 
pect fast exponential decay of any mode, and this system 
is near the boundary of Lyapunov stability. 



i. Does IBGK with overrelaxation collisions approx- 
imate the BGK equation? The BGK equation as well 
as its discrete velocity version (JTJ has a direction of fast 
contraction n*(/) — /. The discrete chain (|3~4")) with (3 
close to 1 has nothing similar. Hence, the approximation 
of a genuine BGK solution by an LBGK chain may be 
possible only if both the BGK and the LBGK chain tra- 
jectories belong to a slow manifold with high accuracy. 
This implies significant restrictions on initial data and 
on the dynamics of the approximated solution, as well as 
fast relaxation of the LBGK chain to the slow manifold. 

The usual Taylor series based arguments from [32] are 
valid for h <C r. If we assume h 3> t, (St ^ A in the 
notation of [32]) then Eqn. (10) of [32] transforms (in 
our notation) into f(x + vh, v,t + h) = f^(x + vh, v,t + 
h) + 0{t) with M = m(f(x + vh,v,t + h)). That is, 
f(x,v,t + h) = flj(x, v,t + h) + 0(t). According to this 
formula, / should almost be at quasiequilibrium after a 
time step /i> t, with some correction terms of order 
t. This first-order in r correction is, of course, the first 
term of the Chapman-Enskog expansion (fTO|) : rfj^) = 
rA/^. (with possible error of order 0(rh)). This is a very 
natural result for an approximation of the BGK solution, 
especially in light of the Chapman-Enskog expansion 0, 
|32| . but it is not the LBM scheme with overrelaxation. 

The standard element in the proof of second-order ac- 
curacy of the BGK equation approximation by an LBGK 
chain uses the estimation of an integral: for time step h 
we obtain from ([T]) the exact identity 

i pt+h 

fi(x + Vih,t + h) = - I (f* mU) (x) ~ Mx,t'))dt\ 

(36) 

where f* rn ^(x) is the quasiequilibrium state that cor- 
responds to the hydrodynamic fields m(f(x,t')). Then 
one could apply the trapezoid rule for integration to the 
right-hand side of (f3"6"l) . The error of the trapezoid rule 
has the order 0(h 3 ): 

f t+h h h 3 ■■ 

J Q(t') dt> = -(Q(t) + Q(t + h)) - -Q(t>), 

where t' S [t,t + h] is a priori unknown point. But for 
the singularly perturbed system |T]), the second deriva- 
tive of the term f* m ^(x) — fi(x,t') on the right hand 

side of (|3"6")l could be of order 1/r 2 , and the whole er- 
ror estimate is 0(/i 3 /t 3 ). This is not small for h > t. 
For backward or forward in time estimates of the inte- 
gral (|36p. errors have the order 0(h 2 /T 2 ). Hence, for 
overrelaxation with h 3> r this reasoning is not appli- 
cable. Many simple examples of quantitative and quali- 
tative errors of this approximation for a singularly per- 
turbed system could be obtained by analysis of a simple 
system of two equations: x — ^(<f)(y) — x), y — ip(x, y) for 
various <fi and tf>. There are examples of slow relaxation 
(instead of fast), of blow-up instead of relaxation or of 
spurious oscillations, etc. 

Hence, one cannot state that LBGK with overrelax- 
ation collisions approximates solutions of the BGK equa- 
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tion. Nevertheless, it can do another job: it can approxi- 
mate solutions of the macroscopic transport equation. As 
demonstrated within this section, the LBGK chain (|34jl . 
after some initial relaxation period, provides a second- 
order approximation to the transport equation, if it goes 
close to the invariant film up to the order 0(h 2 ) (this 
initial relaxation period may have the order 0(h? /t)). 
In other words, it gives the required second-order ap- 
proximation for the macroscopic transport equation un- 
der some stability conditions. 



III. STABILITY AND STABILISATION 
A. Instabilities 

j. Positivity loss. First of all, if / is far from the 
quasiequilibrium, the state Iq (/) (|33|) may be nonphysi- 
cal. The positivity conditions (positivity of probabilities 
or populations) may be violated. For multi- and infinite- 
dimensional problems it is necessary to specify what one 
means by far. In the previous section, / is the whole state 
which includes the states of all sites of the lattice. All the 
involution operators with classical entropies are defined 
for lattice sites independently. Violation of positivity at 
one site makes the whole state nonphysical. Hence, we 
should use here the £oo-norm: close states are close uni- 
formly, at all sites. 

k. Large deviations. The second problem is nonlin- 
earity: for accuracy estimates we always use the as- 
sumption that / is sufficiently close to quasiequilibrium. 
Far from the quasiequilibrium manifold these estimates 
do not work because of nonlinearity (first of all, the 
quasiequilibrium distribution, /jjj, depends nonlinearly 
on M and hence the projection operator, II*, is nonlin- 
ear). Again we need to keep the states not far from the 
quasiequilibrium manifold. 

I. Directional instability. The third problem is a di- 
rectional instability that can affect accuracy: the vector 
/ — IT*(/) can deviate far from the tangent to q (Fig. [TJ. 
Hence, we should not only keep / close to the quasiequi- 
librium, but also guarantee smallness of the angle be- 
tween the direction / — n*(/) and tangent space to q. 

One could rely on the stability of this direction, but 
we fail to prove this in any general case. The directional 
instability changes the structure of dissipation terms: the 
accuracy decreases to the first-order in time step h and 
significant fluctuations of the Prandtl number and viscos- 
ity, etc. may occur. This carries a danger even without 
blow-ups; one could conceivably be relying on nonreliable 
computational results. 

m. Direction of neutral stability. Further, there ex- 
ists a neutral stability of all described approximations 
that causes one-step oscillations: a small shift of / in the 
direction of Ap does not relax back for 3=1, and its 
relaxation is slow for 3 ~ 1 (for small viscosity). This 
effect is demonstrated for a chain of mirror reflections in 
Fig.H 



Invariant film 




FIG. 1: Directional instability: after several iterations the tra- 
jectory is not tangent to the invariant film with the required 
accuracy. 




FIG. 2: Neutral stability and one-step oscillations in a se- 
quence of reflections. Bold dotted line - a perturbed motion, 
A - direction of neutral stability. 

B. Dissipative recipes for stabilisation 

n. Positivity rule. There is a simple recipe for pos- 
itivity preservation [ill . l5o| : to substitute nonpositive 
Iq (f)(x) by the closest nonnegative state that belongs 
to the straight line 

{A/(x) + (l-A)n*(/(x))|AeR} (37) 

defined by the two points, f(x) and its corresponding 
quasiequilibrium state. This operation is to be applied 
point-wise, at the points of the lattice where the positiv- 
ity is violated. The coefficient A depends on x too. Let 
us call this recipe the positivity rule (Fig. [3]) ; it preserves 
positivity of populations and probabilities, but can affect 
the accuracy of approximation. The same rule is neces- 
sary for ELBGK (|35|) when a positive "mirror state" / 
with the same entropy as / does not exists on the straight 
line (|37|) . 

The positivity rule saves the existence of positive so- 
lutions, but affects dissipation because the result of the 
adjusted collision is closer to quasiequilibrium. There 
is a family of methods that modify collisions at some 
points by additional shift in the direction of quasiequi- 
librium. The positivity rule represents the minimal nec- 
essary modification. It is reasonable to always use this 
rule for LBM (as a "salvation rule'' ) . 

o. Ehrenfests' regularisation. To discuss methods 
with additional dissipation, the entropic approach is very 
convenient. Let entropy S(f) be defined for each popu- 
lation vector / = (fi) (below, we use the same letter S 
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FIG. 3: Positivity rule in action. The motions stops at the 
positivity boundary. 

for local-in-space entropy, and hope that the context will 
make this notation clear). We assume that the global 
entropy is a sum of local entropies for all sites. The local 
noncquilibrium entropy is 

AS(f) = S(f*) - S(f), (38) 

where /* is the corresponding local quasiequilibrium at 
the same point. 

The Ehrenfests ' regularisation [ill [l2| provides "en- 
tropy trimming" : we monitor local deviation of / from 
the corresponding quasiequilibrium, and when AS(f, x) 
exceeds a pre-specified threshold value 5, perform local 
Ehrenfests' steps to the corresponding equilibrium. So 
that the Ehrenfests' steps are not allowed to degrade the 
accuracy of LBGK it is pertinent to select the k sites with 
highest AS > 8. The a posteriori estimates of added dis- 
sipation could easily be performed by analysis of entropy 
production in Ehrenfests' steps. Numerical experiments 
show (see, e.g., [ill, EH and Sect. ITy)) that even a small 
number of such steps drastically improves stability. 

To avoid the change of accuracy order "on average" , 
the number of sites with this step should be (D(NS X /L) 
where N is the total number of sites, 6 X is the step of the 
space discretization and L is the macroscopic character- 
istic length. But this rough estimate of accuracy in aver- 
age might be destroyed by concentrations of Ehrenfests' 
steps in the most nonequilibrium areas, for example, in 
boundary layers. In that case, instead of the total num- 
ber of sites N in the estimate Q(NS X /L ) w e should take 
the number of sites in a specific region [59|]. The effects 
of concentration could be easily analysed a posteriori. 

p. Entropic steps for nonentropic equilibria. If the 
approximate discrete equilibrium /* is nonentropic, we 
can use A5k(/) = — <Sk(/) instead of AS(f), where Sk 
is the Kullback entropy. This entropy, 

5 K(/)=-E/' ln (^)' ( 39 ) 

gives the physically reasonable entropic distance from 
equilibrium, if the supposed continuum system has the 
classical BGS entropy. In thermodynamics, the Kull- 
back entropy belongs to the family of Massieu-Planck- 
Kramers functions (canonical or grand canonical poten- 
tials). One can use (|3"9")l in the construction of Ehrenfests' 
regularisation for any choice of discrete equilibrium. 



We have introduced two procedures: the positivity rule 
and Ehrenfests' regularisation. Both improve stability, 
reduce nonequilibrium entropy, and, hence, nonequilib- 
rium fluxes. The proper context for discussion of such 
procedures are the flux-limiters in finite difference and 
finite volume methods. Here we refer to the_ classi- 
cal flux-corrected transport (FCT) algorithm [10( that 
strictly maintains positivity, and to its further develop- 
ments @, in, nn . 

q. Smooth limiters of nonequilibrium entropy. The 
positivity rule and Ehrenfests' regularisation provide 
rare, intense and localised corrections. Of course, it 
is easy and also computationally cheap to organize 
more gentle transformations with smooth shifts of higher 
nonequilibrium states to equilibrium. The following regu- 
larisation transformation distributes its action smoothly: 

/^/* + </,(AS(/))(/-f). (40) 

The choice of the function <fi is highly ambiguous, for ex- 
ample, (j> = 1/(1 + aAS k ) for some a > 0, k > 0. There 
are two significantly different choices: (i) ensemble- 
independent tf> (i.e., the value of <fi depends on the lo- 
cal value of AS only) and (ii) ensemble-dependent </>, for 
example 

1 + (A^/(aE(A^))) fc ~ 1 / 2 
' ~ l + {AS/{aE(AS))) k 

where E(AS') is the average value of AS in the compu- 
tational area, k > 1 and a > 1. It is easy to select 
an ensemble-dependent cf> with control of total additional 
dissipation. 

r. ELBGK collisions as a smooth limiter. On the 
basis on numerical tests, the authors of f5(| claim that 
the positivity rule provides the same results (in the sense 
of stability and absence/presence of spurious oscillations) 
as the ELBGK models, but ELBGK provides better ac- 
curacy. 

For the formal definition of ELBGK (f3"5|) our tests do 
not support claims that ELBGK erases spurious oscil- 
lations (see Sect. IIVI below). Similar observations for 
Burgers equation has been reported in 0]. We under- 
stand this situation in the following way. The entropic 
method consists of at least three components: 

1. entropic quasiequilibrium defined by entropy max- 
imisation; 

2. entropy balanced collisions ([55)) that have to pro- 
vide proper entropy balance; 

3. a method for the solution of the transcendental 
equation S(f) = S(f) to find a = a{f) in ([55]). 

It appears that the first two items do not affect spuri- 
ous oscillations at all, if we solve equation for a(f) with 
high accuracy. Additional viscosity could, potentially, be 
added by use of explicit analytic formulas for a(f). In 
order not to decrease entropy, the errors in these formu- 
las always increase dissipation. This can be interpreted 
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as a hidden transformation of the form (|40| . where the 
coefficients of <j> also depend on /*. 

Compared to flux limitcrs, noncquilibrium entropy lim- 
iters have a great benefit: by summation of all entropy 
changes we can estimate the amount of additional dissi- 
pation the limiters introduce into the system. 



C. Non-dissipative recipes for stabilisation 

s. Microscopic error and macroscopic accuracy. 
The invariant film q is an invariant manifold for the free- 
flight transformation and for the temporal and entropic 
involutions. The linear involution Jo, as well as the EL- 
BGK involution Ie, transforms a point / € q into a point 
/' with / - /' - oA n .(/) + o(J - IT(/)), i.e., the vector 
/ — /' is "almost tangent" to q, and the distance from /' 
to q has the order 0(\\f - n*(/)|| 2 ). 

Hence, if the initial state belongs to q, and the distance 
from quasicquilibrium is small enough (~ 0(h)), then 
during several steps the LBGK chain will remain near q 
with deviation ~ 0(h 2 ). Moreover, because errors pro- 
duced by collisions (deviations from q) have zero macro- 
scopic projection, the corresponding macroscopic error in 
M during several steps will remain of order 0(h 3 ). 

To demonstrate this, suppose the error in /, Sf, is 
of order 0(h k ), and m(Sf) = 0, then for smooth fields 
after a free-flight step an error of higher order appears 
in the macroscopic variables M: m(9ft,(i5/)) = 0(h k+1 ), 
because m{@h{8f)) — m{{@h — !)(£/)) an d ®h — 1 = 
0(h). The last estimate requires smoothness. 

This simple statement is useful for the error analysis 
we perform. We shall call it the lemma of higher macro- 
scopic accuracy: a microscopic error of order 0(h k ) in- 
duces, after a time step h, a macroscopic error of order 
0(h k+1 ), if the field of macroscopic fluxes is sufficiently 
small (here, the microscopic error means the error that 
has zero macroscopic projection). 

t. Restarts and approximation of A f* . The prob- 
lem of nondissipative LBM stabilisation we interpret as 
a problem of appropriate restart from a point that is 
sufficiently close to the invariant film. If h = r and 
collisions return the state to quasiequilibrium, then the 
state belongs to q for all time with high accuracy. For 
h t, formulas for restarting are also available: one can 
choose between and, more flexibly, (f2"3"| . Neverthe- 
less, many questions remain. Firstly, what should one 
take for Af« ? This vector has a straightforward diffcr- 

J M rnn / 

ential definition (110]) (let us also recall that tA ^ * is the 
first Chapman-Enskog noncquilibrium correction to the 
distribution function (|19p ). But numerical differentiation 
could violate the exact-in-space free-flight transformation 
and local collisions. There exists a rather accurate cen- 
tral difference approximation of A f* on the basis of free- 
flight: 

A« f = i(Aj r +A^,) + 0(^) 1 (41) 



where 

A^ =l(Q h (f* M )-Il*(Q h (r M ))), 

&r M = -\{®-h{r M ) Ti*{Q_ h {f* M ))). 

There are no errors of the first-order in (|4"Tj) . The forward 
(At, ) and backward (AT* ) approximations are one or- 

der less accurate. The computation of A.% is of the same 

J M 

computational cost as an LBGK step, hence, if we use the 
restart formula (f22|) with central difference evaluation of 
A/» (|4"Tj) , then the computational cost increases three 
times (approximately) . Non- locality of collisions (restart 
from the distribution f^ (|22|) with a nonlocal expres- 
sion for A f* ) spoils the main LBM idea of exact linear 
free- flight and local collisions: nonlocality is linear and 
exact, nonlinearity is local |53|). One might also consider 
the inclusion of other finite difference representations for 
A/* into explicit LBM schemes. The consequences of 
this combination should be investigated. 

u. Coupled steps with quasiequilibrium ends. The 
mean viscosity lemma allows us to combine different 
starting points in order to obtain the necessary macro- 
scopic equations. From this lemma, it follows that 
the following construction of two coupled steps with 
restart from quasiequilibrium approximates the macro- 
scopic equation ^ with second-order accuracy in time 
step h. 

Let us take f^j as the initial state with given M, then 
evolve the state by 0^, apply the incomplete temporal 
involution/^ ([28]) . again evolve by @h, and finally project 
by n* onto the quasiequilibrium manifold: 

M ^ M' = m(U*(Q h (4(Q h (r M )))))- (42) 

It follows from the restart formula (|22p and the mean 
viscosity lemma that this step gives a second-order in 
time h approximation to the shift in time 2h for (J9j> with 
t = 2(1 - /3)h, 1/2 < J3 < 1. Now, let us replace 1% 
by the much simpler transformation of LBGK collisions 

lo »: 

M i ^ M' = m(IL*(Q h (lP(Q h (r M ))))). (43) 
According to the lemma of higher macroscopic accu- 




FIG. 4: The scheme of coupled steps (35) . 
racy this step (Fig. 2|) also gives a second-order in time 
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h approximation to the shift in time 2h for ([9]) with 
r = 2(1 - 0)h, 1/2 < (3 < 1. The replacement of 7^ 
by 7g introduces an error in / that is of order 0(h 2 ), 
but both transformations conserve the value of macro- 
scopic variables exactly. Hence (due to the lemma of 
higher macroscopic accuracy) the resulting error of cou- 
pled steps in the macroscopic variables M is of order 
0(h 3 ). This means that the method has second-order ac- 
curacy. 

Let us enumerate the macroscopic states in (j43|) : Mo = 
M, Mi/2 = m(<d h (f* M )) and M x = M' . The shift from 
Mo to M 1 / 2 approximates the shift in time h for © with 
t = h. If we would like to model ([9]) with r <C h, then 
t = h means relatively very high viscosity. The step 
from Mi/2 to M\ has to normalize viscosity to the re- 
quested small value (compare to antidiffusion in [9T[To|). 
The antidiffusion problem necessarily appears in most 
CFD approaches to simulation flows with high-Reynolds 
numbers. Another famous example of such a problem 
is the filtering-defiltering problem in large eddy simula- 
tion (LES) [52]. The antidiffusion in the coupled steps 
is produced by physical fluxes (by free-flight) and pre- 
serves positivity. The coupled step is a transformation 
Mo i—* Mi and takes time 2h. The middle point M\/2 is 
an auxiliary state only. 

Let us enumerate the microscopic states in (|43j) : fo = 

fhi A/2 = ®h(fii), ft/2 = n*(A/ 2 )) A/2 = ^o(A/2)' 

ft/2 = fx/2 + (i - P)Chn - A°/2)» fi = ©ft(A/ 2 ) h = 

n*(A7) = f%f , where M x — m(/ 1 ~). Here, in the middle 
of the step, we have 4 points: a free-flight shift of the 
initial state (A/2)' ^ ne corresponding quasiequilibrium 
(f 1/2)1 the mirror image (fi/2) of the point fy 2 with 
respect to the centre A°/ 2 j and the state (ft/2) that is 
the image of /^ 2 after homothety with centre A/2 and 
coefficient 2/3—1. 

For smooth fields, the time shift Oh returns fx 12 to 
the quasiequilibrium manifold with possible error of or- 
der 0(h ). For entropic equilibria, the nonequilibrium 
entropy of the state QhQi/2) is of order 0(h A ). This is 
an entropic estimate of the accuracy of antidiffusion: the 
nonequilibrium entropy of fx/2 could be estimated from 
below as C(M)h 2 , where C(M) > does not depend 
on h. The problem of antidiffusion can be stated as an 
implicit stepping problem: find a point / such that 

m(f) = M, (IT - l)(0 h (/)) - 0. (44) 

This antidiffusion problem is a proper two-point bound- 
ary value problem. In a finite-dimensional space the first 
condition includes N independent equations (where N is 
the number of independent macroscopic variables), the 
second allows N degrees of freedom, because the values 
of the macroscopic variables at that end are not fixed and 
Qh(f) could be any point on the quasiequilibrium mani- 
fold) . Shooting methods for the solution of this problem 
looks quite simple: 



• Method A, 

f n+1 = f n + n*(6 h (/ n )) - B h (f n ), (45) 

• Method B, 

/ n+1 = e_ ft (n*(e h (/ n ))) 

+ /it f -n*(e_ h (n*(e h (/ n )))). 

Method A is: shoot from the previous approximation, 
/„, by Q h , project onto quasiequilibrium, IT (6/, (/„)), 
and then correction of /„ by the final point displacement, 
n*(9fc(/n)) - ©ft(/n)' The value of M does not change, 
because m(Tl*(Q h (f n ))) = m(@ h (f n )). 

Method B is: shoot from the previous approximation, 
/„, by Qh, project onto quasiequilibrium, shoot back- 
wards by Q—h, an( l then correction of M using quasiequi- 
libria (plus the quasiequilibrium with required value of 
M, and minus one with current value of M). 

The initial approximation could be /1/2, and n here is 
the number of iteration. Due to the lemma of higher 
macroscopic accuracy, each iteration (1431) or (|4^|) in- 
creases the order of accuracy (see also the numerical test 
in Sec. [IT]). 

The shooting method A (T4"5|) better meets the main 
LBM idea: each change of macroscopic variable is due to 
a free- flight step (because free-flight in LBM is exact), all 
other operations effect nonequilibrium component of the 
distribution only. The correction of M in the shooting 
method B (T4l)|) violates this requirement. 

The idea that all macroscopic changes are projections 
of free-flight plays, for the proposed LBM antidiffusion, 
the same role as the monotonicity condition for FCT [l(| • 
In particular, free-flight never violates positivity. 

If we a find solution / to the antidiffusion problem with 
M = M x / 2 , then we can take f+ /2 = /° /2 + (1 - /?)(/ - 

A/ 2 ), and Mi = m(Ofe(/ ] y 2 )). But even exact solutions 
of (|4"4")) can cause stability problems: the entropy of / 
could be less than the entropy of A/2> anc ^ blow-up could 
appear. A palliative solution is to perform an entropic 
step: to find a such that S(f®, 2 + a(f — f®, 2 )) — S(f 1 ~, 2 ), 

then use /+ /2 = /° /2 + (1 - (3)a(f - /° /2 ). Even for 
nonentropic equilibria it is possible to use the Kullback 
entropy (|39[) for comparison of distributions with the 
same value of the macroscopic variables. Moreover, the 
quadratic approximation to (f59"| will not violate second- 
order accuracy, and does not require the solution of a 
transcendental equation. 

The viscosity coefficient is proportional to r and sig- 
nificantly depends on the chain construction: for the se- 
quence (|2T)1) we have r = (1 — j3)h//3, and for the sequence 
of steps (|4*5)) r = 2(1 — 0) h. For small 1—0 the later gives 
around two times larger viscosity (and for realisation of 
the same viscosity we must take this into account). 

How can the coupled steps method (|4Uf fail? The 
method collects all the high order errors into dissipation. 
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When the high-order errors accumulated in dissipation 
become compatible with the second-order terms, the ob- 
servable viscosity significantly increases. In our numeri- 
cal tests this catastrophe occurs when the hydrodynamic 
fields change significantly on 2-3 grid steps S x (the char- 
acteristic wave length A ~ 35 x )- The catastrophe point 
is the same for the plain coupled steps (|4"3")l and for the 
methods with iterative corrections or (|4l))) . The ap- 
propriate accuracy requires A > 10S X . On the other hand, 
this method is a good solver for problems with shocks (in 
comparison with standard LBGK and ELBGK) and pro- 
duces shock waves with very narrow fronts and almost 
without Gibbs effect. So, for sufficiently smooth fields 
it should demonstrate second-order accuracy, and in the 
vicinity of steep velocity derivatives it increases viscosity 
and produces artificial dissipation. Hence, this recipe is 
nondissipative in the main order only. 

The main technical task in "stabilisation for accuracy" 
is to keep the system sufficiently close to the invariant 
film. Roughly speaking, we should correct the micro- 
scopic state / in order to keep it close to the invariant film 
or to the tangent straight line {/^f+AA/* |A £ E}, where 
M = m(f). But the general accuracy condition (|23|) gives 
much more freedom: the restart point should return to 
the invariant film in projections on the macroscopic vari- 
ables and their fluxes only. A variant of such regularisa- 
tion was de facto proposed and successfully tested in Q • 
In the simplest realisation of such approaches a problem 
of "ghost" variables [53[ can arise: when we change the 
restart (j2"2"|) to (|2"3"|) . neither moments, nor fluxes change. 
The difference is a "ghost" vector. At the next step, the 
introduced ghost component could affect fluxes, and at 
the following steps the coupling between ghost variables 
and macroscopic moments emerges. Additional relax- 
ation times may be adjusted to suppress these nonhydro- 
dynamic ghost variables (20j . 

v. Compromise between nonequilibrium memory and 
restart rules. Formulas (|22|) and (|23[) prescribe a choice 
of restart states f^j. All memory from previous evolu- 
tion is in the macroscopic state, M, only. There is no 
microscopic (or, alternatively, nonequilibrium) memory. 
Effects of nonequilibrium memory for LBM are not yet 
well studied. For LBGK with overrelaxation, these ef- 
fects increase when (3 approaches 1 because relaxation 
time decreases. We can formulate a hypothesis: ob- 
served sub-grid properties of various LBGK realisations 
and modifications for high-Reynolds number are due to 
nonequilibrium memory effects. 

In order to find a compromise between the restart re- 
quirements (|22p . (|23p and nonequilibrium memory ex- 
istence we can propose to choose directions in concor- 
dance with (|22p . (|23p . where the nonequilibrium entropy 
field (|3"8)) does not change in the restart procedure. If 
after a free-flight step we have a distribution / and find 
a corresponding restart state /„(/) due to a global rule, 
then for each grid point x we can restart from a point 

fmU) + a ( x W™(f) ~ where a ( x ) > is a solu " 

tion of the constant local nonequilibrium entropy equa- 



tionS(f^ f} (x) + a(x)(f^ f} (x)-r mU) (x))) = S(f(x)). 
This family of methods allows a minimal nonequilibrium 
memory - the memory about local cntropic distance from 
quasiequilibrium. 

IV. NUMERICAL EXPERIMENT 
A. Velocities and equilibria 

To conclude this paper we report two numerical ex- 
periments conducted to demonstrate the performance of 
some of the proposed LBM stabilisation recipes from 

Sect, nnj 

We choose velocity sets with entropic equilibria and an 
-ff-theorem in order to compare all methods in a uniform 
setting. 

In ID, we use a lattice with spacing and time step 
h = 1 and a discrete velocity set t>2, V3} := {0, —1, 1} 
so that the model consists of static, left- and right-moving 
populations only. The subscript i denotes population 
(not lattice site number) and fx, fa and denote the 
static, left- and right-moving populations, respectively. 
The entropy is S = —H, with 

H = h log(/i/4) + h log(/a) + h log(/ 3 ), 

(see, e.g., [4l|) and, for this entropy, the local quasiequi- 
librium state /* is available explicitly: 

ft = ^((3 U -l) + 2v/l + 3u 2 ), 
/ 3 * = -|((3i/ + l)-2 V / l + 3^) ! 

where 

n:=y~)fi, u:=-VV/i. 

l i 

In 2D, the realisation of LBGK that we use will employ 
a uniform 9-speed square lattice with discrete velocities 
{Vi\ i = 0, 1, . . . 8}: vo — 0, Vi — (cos((t — l)7r/2), sin((i — 
l)7r/2)) for i = 1,2,3,4, v t = V2(cos((i - 5)f + 
f ),sin((i - 5)f + f)) for i = 5,6,7,8. The numbering 
foi fii ■ ■ ■ j fs are for the static, east, north, west, south, 
northeast, northwest, southwest and southeast-moving 
populations, respectively. As usual, the quasiequilibrium 
state, /*, can be uniquely determined by maximising an 
entropy functional 

stfl—i; 

i 

subject to the constraints of conservation of mass and 
momentum [3j: 

jr-w.n(»-V^i)( ,- Uj ) ■ 

(47) 
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Here, the lattice weights, Wi, are given lattice-specific 
constants: W — 4/9, VFx,2,3,4 = 1/9 an d ^5,6,7,8 = 
1/36. The macroscopic variables are given by the ex- 
pressions 

i — ' n ^ — ' 

i i 

As we are advised in Sect. HID in all of the experiments, 
we implement the positivity rule. 



B. Shock tube 

The ID shock tube for a compressible isothermal fluid 
is a standard benchmark test for hydrodynamic codes. 
Our computational domain will be the interval [0, 1] and 
we discretize this interval with 801 uniformly spaced lat- 
tice sites. We choose the initial density ratio as 1:2 so 
that for x < 400 we set n — 1.0 else we set n = 0.5. 

w. Basic test: LBGK. ELBGK and Coupled steps. 
We will fix the kinematic viscosity of the fluid at v = 10~ 9 
We should take j3 = l/(2u + 1) « 1 - 2u for LBGK and 
ELBGK (with or without the Ehrenfests' regularisation) . 
Whereas, for the coupled step regularisation, we should 
take (3=1 — 1/. 

The governing equations for LBGK are 

f i (x + v i ,t + l) = f*(x,t) + (28-l)(f*(x,t) -fi(x,t)). 

(48) 

For ELBGK (f35|) the governing equations are: 



Mx + v t ,t+l) = (l-(3)f*(x,t)+(3fi(x,t), (49) 

with / = (1 — a)f + af*. As previously mentioned, 
the parameter, a, is chosen to satisfy a constant entropy 
condition. This involves finding the nontrivial root of the 
equation 



S((l-a)f + af*) = S(f). 



(50) 



Inaccuracy in the solution of this equation can introduce 
artificial viscosity. To solve (|5"0)) numerically we employ 
a robust routine based on bisection. The root is solved 
to an accuracy of 10~ 15 and we always ensure that the 
returned value of a does not lead to a numerical entropy 
decrease. We stipulate that if, at some site, no nontriv- 
ial root of (|50p exists we will employ the positivity rule 
instead. 

The governing equations for the coupled step regular- 
isation of LBGK alternates between classic LBGK steps 
and equilibration: 



fi{x + Vi,t + l) 

J*( X) t) + (2(3 



N step odd, 



- fi(x,t)), N stcp even, 
(51) 



where N stC p is the cumulative total number of time steps 
taken in the simulation. For coupled steps, only the re- 
sult of a couple of steps has clear physical meaning: this 



couple transforms f*(x,t) that appears at the beginning 
of an odd step to f*(x,t) that appears at the beginning 
of the next odd step. 
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FIG. 5: Density and velocity profile of the 1:2 isother- 
mal shock tube simulation after 400 time steps using (a) 
LBGK (J48J); (b) ELBGK Jl9j; (e) coupled step regularisa- 
tion (|5ip : In this example, no negative population are pro- 
duced by any of the methods so the positivity rule is redun- 
dant. For ELBGK in this example, (|50p always has a non- 
trivial root. 

As we can see, the choice between the two collision for- 
mulas LBGK g5D or ELBGK @H]) does not affect spuri- 
ous oscillation. But it should be mentioned that the en- 
tropic method consists of not only the collision formula, 
but, what is important, includes the provision of special 
choices of quasiequilibrium that could improve stability 
(see, e.g., [17|)- The coupled steps produce almost no 
spurious oscillations. This seems to be nice, but in such 
cases it is necessary to monitor the amount of artificial 
dissipation and to measure the viscosity provided by the 
method (see below). 

x. Ehrenfests' regularisation. For the realisation of 
the Ehrenfests' regularisation of LBGK, which is in- 
tended to keep states uniformly close to the quasiequi- 
librium manifold, we should monitor nonequilibrium en- 
tropy AS ([55)) at every lattice site throughout the sim- 
ulation. If a pre-specified threshold value S is exceeded, 
then an Ehrenfests' step is taken at the corresponding 
site. Now, the governing equations become: 



fi(x + Vi,t + l) 

'f*(x,t) + (2/3- l)(f*(x,t) 
^f*(x,t), otherwise, 



fi(x,t)), AS<6, 



(52) 
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Furthermore, so that the Ehrenfests' steps are not al- 
lowed to degrade the accuracy of LBGK it is pertinent to 
select the k sites with highest AS > 8. The a posteriori 
estimates of added dissipation could easily be performed 
by analysis of entropy production in Ehrenfests' steps. 




FIG. 6: Density and velocity profile of the 1:2 isothermal 
shock tube simulation after 400 time steps using Ehrenfests' 
regularisation (g2j with (a) (k,S) = (4, 10~ 3 ); (b) (k,S) = 
(4, 10~ 4 ). Sites where Ehrenfests' steps are employed are in- 
dicated by crosses. Compare to Fig. [5^. 



1 ■ 

u0.8 
0.6 
0.4 



a) 







1 ■ 

n0.8 
0.6 



b) 



1 

u0.8 
0.6 
0.4 



O 







0.5 

X 



0.5 

X 



0.5 

X 



ES 



1 °0 



ES 



ES 



200 

t 



400 



200 

t 



400 



200 

t 



400 



FIG. 8: LBGK (gSJ regularised with Ehrenfests' steps J52). 
Density profile of the 1:2 isothermal shock tube simulation 
and Ehrenfests' steps histogram after 400 time steps using 
the tolerances (a) (k,S) = (1,10^ 4 ); (b) (k,8) = (4, 10" 4 ); (c) 
(k,5) = (8, 10 -4 ). Sites where Ehrenfests' steps are employed 
are indicated by crosses. Compare to Fig. [5jt. 



1 - 

n0.8 
0.6 

0.4- 



a) 







0.5 

X 



11 
nO.i 
0.6 
0.4 



6 

ES 4 
2 

: °0 
15 

ES™ 
5 



b) 



o 



0.5 

X 



1 



1 ■ 

n0.8 
0.6 
0.4- 



c) 



ES 







0.5 

X 



40 
20 








200 

t 



400 



200 

t 



400 



200 

t 



400 



FIG. 7: LBGK (f48j) regularised with Ehrenfests' steps (|52j) . 
Density profile of the 1:2 isothermal shock tube simulation 
and Ehrenfests' steps histogram after 400 time steps using 
the tolerances (a) (k,S) = (oo, 10" 3 ); (b) (k,6) = (oo, 10 -4 ); 
(c) (k,8) — (oo, 10 -5 ). Sites where Ehrenfests' steps are em- 
ployed are indicated by crosses. Compare to Fig. [S^i. 

In the example in Fig. O we have considered fixed tol- 
erances of (k, 8) = (4, 1CT 3 ) and (fc, 8) = (4, 10" 4 ) only. 



We reiterate that it is important for Ehrenfests' steps 
to be employed at only a small share of sites. To il- 
lustrate, in Fig. [7] we have allowed k to be unbounded 
and let 8 vary. As 5 decreases, the number of Ehrenfests' 
steps quickly begins to grow (as shown in the accompany- 
ing histograms) and excessive and unnecessary smooth- 
ing is observed at the shock. The second-order accuracy 
of LBGK is corrupted. In Fig. EJ we have kept 5 fixed at 
8 = 10~ 4 and instead let k vary. We observe that even 
small values of k (e.g., k = 1) dramatically improves the 
stability of LBGK. 



C. Accuracy of coupled steps 



Coupled steps (|43|) give the simplest second-order ac- 
curate stabilization of LBGK. Stabilization is guaranteed 
by collection of all errors into dissipative terms. But this 
monotone collection of errors could increase the higher 
order terms in viscosity. Hence, it seems to be necessary 
to analyze not only order of errors, but their values too. 

For accuracy analysis of coupled steps we are inter- 
ested in the error in the antidiffusion step (|44|) . We 
analyse one coupled step for f3 = 1. The motion 
starts from a quasiequilibrium /o = /L, then a free- 
flight step fy 2 = Qhiflt), after that a simple reflec- 
tion f 1 / 2 = Io(fi/ 2 ) with respect to the quasiequilib- 
rium centre f®, 2 = U*(f~, 2 ), again a free-flight step, 
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FIG. 9: (a) The £2 estimate of middle point displace- 
ment l(53)) : for coupled steps (|43|) (diamonds), one (trian- 
gles), two (squares) and three (dots) shooting iterations (135)) ; 
(b) The £2 estimate of nonequilibrity of the final point (|54[l : 
for coupled steps (|43[) (diamonds) , one (triangles) , two (stars) 
and three (squares) shooting iterations (|46|) . 



fi = &k(f 1/2)1 an d finally a projection onto quasiequi- 
librium, h = n*(/r). 

In the first of two accuracy tests, two types of errors 
are to be studied. The middle point displacement is 



,5 CS = ||/° /2 - n*(e_ fc (/r))||/||/o - A° /2 || 



(53) 



To estimate nonequilibrity of the final point (i.e., ad- 
ditional dissipation introduced by last projection in the 
coupled step) we should compare the difference /£" — /1 
to the difference at the middle point fy 2 
introduce 



fy 2 . Let us 



ii/r-/iii 2 /ii/r /2 -A /2 i 



(54) 



In our tests (Fig. [9]) we use the £2 norm. 

We take the ID 3-velocity model with entropic equilib- 
ria. Our computational domain will be the interval [0, 1] 
which we discretize with 1001 uniformly spaced lattice 
sites. The initial condition is n(x, 0) = 1 + 0.2 sin (2ttujx), 
u{x, 0) = 0.1 cos (2iru>x) and we employ periodic bound- 
ary conditions. We compute a single coupled step for 
frequencies in the range id = 1, 2, . . . , 1000 (Fig. [9]). 

The solution / to the antidiffusion problem could be 
corrected by the shooting iterations ([45]) and (|46)l . The 
corresponding errors for method B (|46|) are also presented 
in Fig. [91 We use (5 cs ,i and er C sa, for i = 1, 2, . . . , to denote 
each subsequent shooting of (fSU)) and ([5"4"[) , respectively. 



We observe that the nonequilibrity estimate, a cs , 
blows- up around the wavelength l/to ~ 38 x . Simultane- 
ously, the middle point displacement S cs has value around 
unity at the same point. We do not plot the results for 
larger values of uj as the simulation has become mean- 
ingless and numerical aliasing will now decrease these 
errors. The same critical point is observed for each sub- 
sequent shooting as well. For this problem, the shooting 
procedure is demonstrated to be effective for wavelengths 
l/u) < 10S X . 

For the second accuracy test we propose a simple test 
to measure the observable viscosity of a coupled step 
(and LBGK) simulation. We take the 2D isothermal 
9-velocity model with entropic equilibria. Our compu- 
tational domain will a square which we discretize with 
[L + 1) x (L + 1) uniformly spaced points and periodic 
boundary conditions. The initial condition is n(cc, y) = 1, 
ui(x,y) = and U2(x,y) = uq siri(2irx / L) , with uq — 
0.05. The exact velocity solution to this problem is an 
exponential decay of the initial condition: u\(x, y, t) = 0, 
U2{x,y,t) — uoexjp(— Auot /(ReL))sm(2irx/L) , where A 
is some constant and Re = Re((3) = uqL/v{0) is the 
Reynolds number of the flow. Here, v = v((3) is the the- 
oretical viscosity of the fluid: v = 1 — /? for the coupled 
steps gS]) and v = (1//3- l)/2 for LBGK. 

Now, we simulate the flow over L/vq time steps and 
measure the constant A from the numerical solution. We 
do this for both LBGK and the coupled steps (|43[) for 
L = 100 and for L = 200. The results (Fig. [TO]) show us 
that for coupled steps (and for LBGK to a much lesser 
extent) the observed viscosity is higher than the theoret- 
ical estimate, hence the observed Re is lower than the 
estimate. In particular, the lower-resolution (L = 100) 
coupled steps simulation diverges from LBGK at around 
Re = 500. The two times higher-resolution (L = 200) 
simulations are close to around Re = 0(1000), after 
which there begins to be a considerable increase in the 
observable viscosity (as explained within Sect. IIII C)) . 



D. Flow around a square-cylinder 

The unsteady flow around a square-cylinder has been 
widely experimentally investigated in the literature (see, 
e.g., [lj| E(|[57j])- The computational set up for the flow 
is as follows. A square-cylinder of side length L, initially 
at rest, is immersed in a constant flow in a rectangular 
channel of length 30L and height 25L. The cylinder is 
place on the centre line in the y-direction resulting in a 
blockage ratio of 4%. The centre of the cylinder is placed 
at a distance 10. 5L from the inlet. The free-stream ve- 
locity is fixed at (itoo^oo) = (0.05,0) (in lattice units) 
for all simulations. 

On the north and south channel walls a free-slip bound- 
ary condition is imposed (see, e.g., [53]). At the inlet, 
the inward pointing velocities are replaced with their 
quasiequilibrium values corresponding to the free-stream 
velocity. At the outlet, the inward pointing velocities 
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FIG. 10: Numerically computed value of A versus Reynolds 
number for 2D accuracy test: LBGK with L = 100 (dia- 
monds); coupled steps (|43[) with £ = 100 (triangles); LBGK 
with L = 200 (stars); coupled steps (g3| with L = 200 
(squares). 



are replaced with their associated quasiequilibrium values 
corresponding to the velocity and density of the penulti- 
mate row of the lattice. 

y. Maxwell boundary condition. The boundary con- 
dition on the cylinder that we prefer is the diffusive 
Maxwell boundary condition (see, e.g., [lj|), which was 
first applied to LBM in [2| • The essence of the condition 
is that populations reaching a boundary are reflected, 
proportional to equilibrium, such that mass-balance (in 
the bulk) and detail-balance are achieved. We will de- 
scribe two possible realisations of the boundary condi- 
tion - time-delayed and instantaneous reflection of equili- 
brated populations. In both instances, immediately prior 
to the advection of populations, only those populations 
pointing in to the fluid at a boundary site are updated. 
Boundary sites do not undergo the collisional step that 
the bulk of the sites are subjected to. 

To illustrate, consider the situation of a wall, aligned 
with the lattice, moving with velocity M W aii and with out- 
ward pointing normal to the wall pointing in the positive 
y-direction (this is the situation on the north wall of the 
square-cylinder with u wa ii = 0). The time-delayed reflec- 
tion implementation of the diffusive Maxwell boundary 
condition at a boundary site (x, y) on this wall consists 
of the update 



fi(x,y,t + l) = a/*(w wa n), 



2,5,6, 



with 



= h(x,y,t) + f 7 (x,y,t) + f s (x,y,t) 

/|(ltwau) + ft ("wall) + / 6 * ("wall) 

Whereas for the instantaneous reflection implementation 
we should use for a: 

f 4 (x,y + l,t)+ f 7 (x + l,y + l,t)+ f s {x - l,y + l,t) 

/| ("wall) + /| ("wall) + fi («wall) 



Observe that, because density is a linear factor of the 
equilibria (j4"T|) , the density of the wall is inconsequential 
in the boundary condition and can therefore be taken as 
unity for convenience. 

We point out that, although both realisations agree 
in the continuum limit, the time-delayed implementation 
does not accomplish mass-balance. Therefore, instanta- 
neous reflection is preferred and will be the implementa- 
tion that we employ in the present example. 

Finally, it is instructive to illustrate the situation for 
a boundary site {x, y) on a corner of the square-cylinder, 
say the north-west corner. The (instantaneous reflection) 
update is then 



fi(x,y,t+ 1) = /3/*(u W aii), 



i = 2,3, 5, 6, 7, 



where 



P — A) /Avail ; 

Pa = h(x - 1,3/,*) + h(x,y + l,t) 

+f 5 (x - 1, y - 1, t) + } 7 {x + 1, y + 1, t) 
+f 8 (x-l,y + l,t); 

Aral] = f% (Wwall) + fs («wall) 

+f£(Uwau) + /e(«mJ]) + //(""wall)- 

z. Strouhal-Reynolds relationship. As a test of the 
Ehrenfests' regularisation (15"2")) , a series of simulations, 
all with characteristic length fixed at L = 20, were con- 
ducted over a range of Reynolds numbers Re = Lu^jv. 
The parameter pair (k, 6), which control the Ehrenfests' 
steps tolerances, are fixed at (L/2, 10~ 3 ). 

We are interested in computing the Strouhal-Reynolds 
relationship. The Strouhal number St is a dimensionless 
measure of the vortex shedding frequency in the wake of 
one side of the cylinder: St = Lf^/uoa, where f u is the 
shedding frequency. 

For our computational set up, the vortex shedding fre- 
quency is computed using the following algorithmic tech- 
nique. Firstly, the s-component of velocity is recorded 
during the simulation over i max = 1250-L/itoo time 
steps. The monitoring points is positioned at coordi- 
nates (4L, — 2L) (assuming the origin is at the centre of 
the cylinder). Next, the dominant frequency is extracted 
from the final 25% of the signal using the discrete Fourier 
transform. The monitoring point is purposefully placed 
sufficiently downstream and away from the centre line 
so that only the influence of one side of the cylinder is 
recorded. 

The computed Strouhal-Reynolds relationship using 
the Ehrenfests' regularisation of LBGK is shown in 
Fig. [TTJ The simulation compares well with Okajima's 
data from wind tunnel and water tank experiment [46| . 
The present simulation extends previous LBM studies of 
this problem p], 0] which have been able to quantitively 
captured the relationship up to Re = 0(1000). Fig. [TT1 
also shows the ELBGK simulation results from Fur- 
thermore, the computational domain was fixed for all 
the present computations, with the smallest value of 
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FIG. 11: Variation of Strouhal number St as a function of 
Reynolds. Dots are Okajima's experimental data [46l ] (the 
data has been digitally extracted from the original paper). 
Diamonds are the Ehrenfests' regularisation of LBGK and 
the squares are the ELBGK simulation from [l[. 



the kinematic viscosity attained being v = 5 x 1CP 5 at 
Re = 20000. It is worth mentioning that, for this charac- 
teristic length, LBGK exhibits numerical divergence at 
around Re = 1000. We estimate that, for the present 
set up, the computational domain would require at least 
O(10 7 ) lattice sites for the kinematic viscosity to be large 
enough for LBGK to converge at Re — 20000. This is 
compared with O(10 5 ) sites for the present simulation. 



V. CONCLUSIONS 

In this paper, we have analysed LBM as a discrete dy- 
namical system generated in distribution space by free- 
flight for time 5 t = h and involution (temporal or en- 
tropic, or just a standard LBGK reflection that approx- 
imates these involutions with second-order accuracy). 
Dissipation is produced by superposition of this involu- 
tion with a homothety with centre in quasiequilibrium 
and coefficient 2(3 — 1. 

Trajectories of this discrete dynamical system are pro- 
jected on to the space of macroscopic variables, hydro- 
dynamic fields, for example. The projection of a time 
step of the LBM dynamics in distribution space approxi- 
mates a time shift for a macroscopic transport equation. 
We represent the general form of this equation ©, and 
provide necessary and sufficient conditions for this ap- 
proximation to be of second-order accuracy in the time 
step h (j22j) . (|23|) . This analysis includes conditions on 
the free-flight initial state, and does not depend on the 
particular collision model. 

It is necessary to stress that for free-flight the space 
discretization is exact (introduces no errors), if the set of 
velocities consists of automorphisms of the grid. 

It seems natural to discuss the LBM discrete dynamical 
system as an approximate solution to the kinetic equa- 



tion, for example, to the BGK kinetics with a discrete 
velocity set (JT|). With this kinetic equation we introduce 
one more time scale, r. For h > t (overrelaxation) the 
discrete LBM does not give a second-order in time step 
h approximation to the continuous-in-time equation {1} . 
This is obvious by comparison of "fast" direction relax- 
ation times: it is r for |T]) and h/(2(l — /?)) ~ h 2 /r 
for discrete dynamics (see also [58]). Nevertheless, the 
"macroscopic shadow" of the discrete LBM with overre- 
laxation approximates the macroscopic transport equa- 
tion with second-order in time step h accuracy under the 
conditions (J22|) and (f2"3"]) . 

We have presented the main mechanisms of observed 
LBM instabilities: 

1. positivity loss due to high local deviation from 
quasiequilibrium ; 

2. appearance of neutral stability in some directions 
in the zero viscosity limit; 

3. directional instability. 

We have found three methods of stability preservation. 
Two of them, the positivity rule and the Ehrenfests' reg- 
ularisation, are "salvation" (or "SOS") operations. They 
preserve the system from positivity loss or from the lo- 
cal blow-ups, but introduce artificial dissipation and it 
is necessary to control the number of sites where these 
steps are applied. In order to preserve the second-order 
of LBM accuracy, in average, at least, it is worthwhile 
to perform these steps on only a small number of sites; 
the number of sites should not be higher than 0{N5 X / L), 
where ./V is the total number of sites, L is the macroscopic 
characteristic length and S x is the lattice step. Moreover, 
because these steps have a tendency to concentrate in the 
most nonequilibrium regions (boundary layers, shock lay- 
ers, etc.), instead of the total number of sites one can use 
an estimate of the number of sites in this region. 

The positivity rule and the Ehrenfests' regularisation 
are members of a wide family of "nonequilibrium entropy 
limiters" that will play the same role, for LBM, as the 
flux limiters play for finite difference, finite volume and 
finite element methods. We have described this family 
and explained how to use entropy estimates for nonen- 
tropic equilibria. The great benefit of the LBM methods 
is that the dissipation added by limiters could easily be 
estimated a posteriori by summarising the entropy pro- 
duction. 

Some practical recommendation for use of nonequilib- 
rium entropy limiters are as follows: 

• there exists a huge freedom in the construction of 
these limiters; 

• for any important class of problems a specific opti- 
mal limiter could be found; 

• one of the simplest and computationally cheapest 
nonequilibrium entropy limiters is the Ehrenfests' 
regularisation with equilibration at k sites with 
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highest nonequilibrium entropy AS > S (the (k, 6)- 
rule); 

• the positivity rule should always be implemented. 

The developed restart methods (Sec. IIII CI) (including 
coupled steps with quasicquilibrium ends) could provide 
second-order accuracy, but destroy the memory of LBM. 
This memory emerges in LBM with overrelaxation be- 
cause of slow relaxation of nonequilibrium degrees of free- 
dom (there is no such memory in the continuous-in-time 
kinetic equation with fast relaxation to the invariant slow 
Chapman-Enskog manifold) . Now, we have no theory of 
this memory but can suggest a hypothesis that this mem- 
ory is responsible for the LBM sub-grid properties. A 
compromise between memory and stability is proposed: 
one can use the directions of restart to precondition col- 
lisions, and keep the memory in the value of the field of 
local nonequilibrium entropy AS (or, for systems with 
nonentropic equilibria, in the value of the corresponding 
Kullback entropies ([39]) ). Formally, this preconditioning 
generates a matrix collision model [53| with a specific 
choice of matrix: in these models, the collision matrix is 
a superposition of projection (preconditioner), involution 
and homothety. A return from the simplest LBGK colli- 
sion to matrix models has been intensively discussed re- 
cently in development of the multirelaxation time (MRT) 
models (for example, [2l|, [3||, see also [13] for matrix 
models for modelling of nonisotropic advection-diffusion 
problems, and [44j for regularisation matrix models for 
stabilisation at high- Reynolds numbers). 

For second-order methods with overrelaxation, ade- 
quate second-order boundary conditions have to be de- 



veloped. Without such conditions either additional dis- 
sipation or instabilities appear in boundary layers. The 
proposed schemes should now be put through the whole 
family of tests in order to find their place in the family 
of the LBM methods. 

Recently, several approaches to stable LBM modelling 
of high- Reynolds number flows on coarse grids have been 
reported [2l|, [2^, [38| . Now it is necessary to understand 
better the mechanisms of the LBM sub-grid properties, 
and to create the theory that allows us to prove the ac- 
curacy of LBM for under-resolved turbulence modelling. 
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